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THEORETICAL MODELING OF PROPAGATION OF MAGNETO-ACOUSTIC WAVES IN MAGNETIC 

REGIONS BELOW SUNSPOTS. 
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ABSTRACT 

We use 2D numerical simulations and eikonal approximation, to study properties of MHD waves 
traveling below the solar surface through the magnetic structure of sunspots. We consider a series 
of magnetostatic models of sunspots of different magnetic field strengths, from 10 Mm below the 
photosphere to the low chromosphere. The purpose of these studies is to quantify the effect of the 
magnetic field on local helioseismology measurements by modeling waves excited by sub-photospheric 
sources. Time-distance propagation diagrams and wave travel times are calculated for models of 
various field strength and compared to the non- magnetic case. The results clearly indicate that the 
observed time-distance helioseismology signals in sunspot regions correspond to fast MHD waves. The 
slow MHD waves form a distinctly different pattern in the time-distance diagram, which has not been 
detected in observations. The numerical results are in good agreement with the solution in the short- 
wavelength (eikonal) approximation, providing its validation. The frequency dependence of the travel 
times is in a good qualitative agreement with observations. 

Subject headings: MHD; Sun: magnetic fields; Sun: oscillations; Sun: helioseismology 



1. INTRODUCTION 



Local helioseismology (time-distance helioseismology, 
acoustic holography and ring-diagram analysis) provides 
valuable information about the physical properties of so- 
lar sub-surface lavers fe.g.lDuy all et al. 1993; K osovichev 
19991. 120021: iKosovichev et al.ri2000: Zhao Kosovichev 



20031: iBraun fc LindsevI l2000l : [Hill 1 19881 : iHaber et al 



2000 ). Time-distance he lioseismology (|Duvall et al 



19931 : IKosovichev fc Duvalll [T997) makes use of wave 
travel times measured for wave packets traveling between 
various points on the solar surface through the interior. 
By inversion of these measurements, variations of the 
wave speed and velocities of mass flows can be recovered 
below the visible solar surface. Inversion results for the 
travel times have been obtained for quiet Sun regions as 
well as for magnetic active regions including sunspots. In 
most of these cases, the variations of the travel times were 
assumed to be caused by mass flows and magneto-sonic 
speed fluctuations along the wave paths. It is known 
that sunspots possess strong magnetic fields with a com- 
plicated structure in the visible layers of the Sun where 
the Doppler measurements used by local helioseismol- 
ogy are taken (|Solankill2003l ). Consequently, such mag- 
netic field can cause important effects on helioseismic 
waves, beyond the perturbation theo ries emp loyed for 
helioseismic data analysi s (Kosovichev fc Duva h 1997; 
iBirch fc KosovichevI l2000l : iGizon fc BirchI 120021 ) . Mag- 
netic field can change the acoustic excitation rate and 
spectral properties of solar oscillations (see the recent 
work bv' Jacoutot et al.ll200 8'): produce new wave modes; 
change wave propagation speeds; change the wave prop- 
agation paths and reflection characteristics at the near- 
surface layers. 
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The influence of magnetic field on the interpreta- 
tion of local helioseismology measurements has not 
been full y explored. T h eoreti ca l effor t s hav e been 
made bv I Crouch fc Callyl (I2003D. [Caliyl (|2QQ^ , 
Schunker fc Callyl (I2006D . ICallv fc GoossensI 



M) 



mm 



Schunker et al.l ( 20081 ) to include mode conversion and 



model waves in magnetized structures by means of 
analytical and ray-path theories. These studies con- 
firm the potential importance of magnetic effects. A 
more complete understanding of the problem can be 
reached by direct forward modeling of helioseismic 
data, by solving numerically the MHD equation s in 



20061: iKhomen 


ko fc ColladosI l2006l: Hanaso^e 


20081: 


Cameron et al.l 


20071: iParchevskv fc Kosovichev 


|2008|). 



In these papers, the MHD equations, governing the prob- 
lem, are solved for magnetic field configurations resem- 
bling a sunspot or a part of a sunspot. The results 
have demonstrated how magnetic fields change the speed, 
amplitude and spectral properties of the waves, and 
also produce wave transformation and refraction in the 
sunspot atmosphere. 

Observational helioseismology over active regions 
has demonstrated repeatedly that waves travel faster 
through su nspots in compariso n with quiet Sun measure- 
ments (e.g. iDuvah et al.lll993l : iBr audi 199 7f ). This fact is 
still the subject of present investigations and its explana- 
tion is still far from clear. Although, given that sunspots 
possess an intense field at photospheric levels, direct or 
indirect the influence of the magne tic field is generally ac- 
cepted (e.g. iHindman et aDll997l ). At sub-photospheric 
levels the influence of the magnetic field becomes less 
and less important with depth as pressure forces domi- 
nate over magnetic forces. This is why the phenomena 
related to magnetic field are referred to as "surface ef- 
fects". Explanations base d on wave perturbations due 
to scattering phase shifts ( Braunin"997l : iLindsey fc BraunI 
[2005a. b) taking place within a few Mm below the pho- 
tosphere or local wave speed variations below sunspots 
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(|KQsovichev et al.ll2QQQl ) have been suggested to explain 
the different travel times measured in sunspots and in 
the quiet photosphere. 

In this paper, we study wave properties and variations 
of travel times with the help of numerical simulations 
and Eikonal solutions for wave propagation in sunspot- 
like magnetic structures in stratified surface and sub- 
surface layers. We perform 2D numerical simulations 
of the waves pr oduced by a single source for a series of 
model sunspots ( Khomenko fc Col lados 2008). We study 
various types of MHD waves and their propagation prop- 
erties numerically. In addition, for the fast magneto- 
acoustic mode we obtain a solution in the eikonal ap- 
proximation, which is often impor tant for cor rect inter- 
pretation of helioseismic results (Gough 2007). Our nu- 
merical results clearly indicate that the fast mode is the 
primary source of helioseismic signals. We study the fre- 
quency dependence of the wave travel times and com- 
pare the results of the eikonal approach with our nu- 
merical simulations and observational data. Our results 
provide insight into the properties of helioseismic waves 
in magnetic regions and are important for interpretation 
of time-distance helioseismology measurements. 

2. NUMERICAL PROCEDURE 

We us e the numerical simulation code described previ- 
ously by iKhomenko fc ColladosI (|2006l ): lKhomenko et al.l 
(|2QQ8l ). The code solves the basic non-linear equations 
of the ideal MHD, written in conservative form: 



dp 

dt 



d{pV) 
dt 



V ■ [pVV + (P - 
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(2) 



— + y.[{E+P+ — )V-B{^)]=pV-g+pQ, (3) 



(4) 



where I is the identity tensor and E is the total energy: 



7- 



1 



Stt 



(5) 



All other symbols have their usual meaning. Gravity has 
been taken constant, with a value of 2.74x10^ cm/sec^. 
The term describing radiative energy losses, Q, is set to 
zero in the present study. A Perfectly Matc hing; Layer 
(PML) boundary condition (jBerengerl 11994 ) is applied 
to all boundaries of the simulation domain with the aim 
of minimizing wave reflection effects at the boundaries. 
This performs rather well in our simulations. The thick- 
ness of the PML layer in our simulations is 20 grid points 
(1 Mm). 

The code solves the non-linear MHD equations for per- 
turbations, that are obtained by subtracting the equa- 
tions of initial magnetostatic equilibrium from Eqs. [TJHl 
However, for the purpose of the present study, we keep 
the perturbation amplitudes as low as possible, making 
the 2nd order non-linear terms relatively unimportant. 
In this paper, we use the 2D version of the code. 

The wave source is modeled as a time-dependent 
term S{x^ z, t) in the right hand side of the momentum 



equation (Eg. [2]) in the same w ay as in the paper by 
iParchevskv fc KosovichevI (|2008r ): 



f{x,Z,T) = 




n 2 



^src 



(1 - 2r^) e-^' if r < Rsrc 

if r > Rsrc, 
(6) 

where Rsrc is the source radius, taken to be 250 
km; r is the distance from the source center, r = 
^/{x — Xsrc^ + — ^src)^- The parameter r is given by 
the equation 

uj(t — to) 47r 
r= \ to<t<to + — , (7) 

Z UJ 

where co is the central, dominating, frequency of the 
source and to is the source start moment. In our model 
the source represents the vertical force in the momen- 
tum equation, S = (0, f)^. The central frequency of the 
source is 3.3 mHz. In all the simulations described here, 
the center of the source is located at depth —200 km, 
below the photospheric level of the reference model S of 
IChristensen-Dalsgaard et al] (|l996f ). The depth location 
of the source in our simulations corresponds to the recent 
studies of the excita tion of waves in magnetic regions by 
IJacoutot et al.' ^2008). 

The source generates a broad spectrum of acoustic and 
surface gravity waves in the 5-minute range, resembling 
the wave excitation in the real Sun. In the case of a 
sunspot model, these modes are modified by the presence 
of magnetic field. 

3. SUNSPOT MODEL AND SIMULATION SETUP 

The sunsp ot models are calculated using the method 
proposed in IKhomenko fc ColladosI (|2008l ). Here we 
use three magnetostatic models of sunspots with simi- 
lar magnetic field structure but different values of the 
magnetic field strength. In this way we are able to study 
the effects of the magnetic field strength on the propa- 
gating waves apart from the effects of the magnetic field 
inclination. Some characteristics of these models are dis- 
played in Figures [T] and O The magneto-static equilib- 
rium is calculated in the cylindrical coordinates assuming 
azimuthal asymmetry and no twist of the magnetic field 
lines. At large distances from the axis, the sunspot mod- 
els merge smoothly into the quiet Sun model, fo r which 
we t ake the model S of Christensen-Dalsga ard et al.l 
(|l996l ). modified to make it stab le with respect to c onvec- 
tive motions (see P archevskv fc KosovichevI l2007l ). The 
sound speed variations along the radial distance are ev- 
ident from Fig. O The sunspot models have a low- 
temperature zone below the surface down to, about, —2 
Mm depth. Below this depth the temperature gradient in 
the horizontal direction is small and no high-temperature 
zone is introduced in the present study. Note that the 
ratio between the acoustic and the Alfven speeds squared 
changes by orders of magnitude with height, and that the 
gradient of this ratio is largest in the surface layers. 

For the wave simulations, the models are transformed 
from cylindrical coordinates (r, 0, z) to 3D cartesian co- 
ordinates. In a general case, all three components of the 
magnetic field vector are non-zero in the cartesian co- 
ordinates. In the case of the 2D simulations, we chose 
a slice where the By component of the field is equal to 
zero and calculate the wave propagation in this plane. 
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Fig. 1. — Top panel: topology of the magnetic field lines for 
the three sunspot models with -Bp^o^ 

=0.9 kG (red lines), 1.5 kG 
(green lines) and 2.4 kG (blue lines). The horizontal solid lines 
mark the levels where cs = va- The contours of the magnetic 
field strength of B = 1 kG are shown by dotted curves for each 
case. Two yellow dots mark the locations of the wave source in the 
two sets of our simulations. Middle panel: enlarged view of the 
upper part of the sunspot models marked by the box in the top 
panel. The vertical axis has been expanded for better visualization. 
Dashed lines are the contours of constant optical depth logrs = 
— 1.6. Bottom panels: the height dependence of the magnetic field 
strength and the gas pressure at the axis for these models. 
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Fig. 2. — The distribution with the radial distance (left panels) 
and with the depth (right panels) of the magnetic field strength, 
B] sound speed, c^; Alfven speed, va] their ratio, Cg/v^; and the 
acoustic cut-off frequency, ljc (calculated in the assumption of an 
isothermal atmosphere) for the Bphot=2.4 kG sunspot. 



By adopting this approach we exclude Alfven waves from 
consideration. 

The simulation domain is 40 Mm in the horizontal 
direction and 10.5 Mm in the vertical direction. The 
sunspot axis is located at the center of the simulation 
domain, and the upper 0.5 Mm of the domain are lo- 
cated above Z = km height. Above Z = 0.5 Mm we 
introduce an isothermal PML layer, 1 Mm thick. The 
grid resolution of the simulations is 50 km in the vertical 
and 100 km in the horizontal directions. The temporal 
resolution depends on the CFL numerical stability con- 
dition, limited by the rapid growth of the Alfven speed, 
Va, with height, changing typically between 0.01 and 0.08 
seconds. For each of the three sunspot models we per- 
form two simulation runs with the wave source located 
at Xq = 6 Mm and 12 Mm to the left of the sunspot axis. 

In addition to the simulations of waves in the sunspot 
models, we have also performed a reference simulation for 
a model without magnetic field and without horizontal 
variations in thermodynamic parameters. The vertical 
stratification of the thermodynamic parameters in this 
case corresponds to the modified model S. 

Several recent works on helioseismology point out the 
presence of the so-called "surface effe cts" in the mea- 
surements in active regions (se e e.g. iL indsev fc BraunI 
2005a,b; Schunker et al. 2005; Zhao fc K osovichev 2006; 
[Rajaguru et al.l I2007D . Due to the Wilson depression 
in sunspots the formation heights of spectral lines are 
altered introducing additional phase shifts into the ve- 
locity measurements, as compared to non-magnetic re- 
gions. Our sunspot models describe the Wilson depres- 
sion effect since the whole atmosphere inside the mag- 
netic concentration had been shifted downwards. In or- 
der to take these "surface effects" into account, we have 
computed the optical depth scale in our sunspot mod- 
els. Th is has been done with the help of SIR code 
(|Ruiz Cobo fc del Toro Iniestalll992[ ). As expected, the 
Ts scale varies with radial distance from sunspot center. 
This variation can be appreciated at the middle panel of 
Fig. [H where the contours of the constant optical depth 
logTs = — 1.6 are given. The optical depth logrs = — 1.6 
represents the layer where typical photospheric spectral 
lines are formed. In the rest of the paper we call "photo- 
spheric velocity" the velocity measured at logrs = — 1.6 
level. Note that the height variations in the position of 
this layer can be as large as 0.2-0.4 Mm, depending on 
the distance from the axis in the sunspot models. In the 
non-magnetic simulation the "photospheric velocity" at 
logTs = — 1.6 always corresponds to the same height, as 
there is no Wilson depression. 

4. WAVE FIELD FROM A SINGLE SOURCE 

In this section we describe the results of the simula- 
tions for the three model sunspots and the quiet Sun 
model for two positions of the oscillation source. In the 
first series of simulations, the source is located in the 
"umbral" part of the model at a distance of 6 Mm from 
the sunspot axis, where the field is relatively strong. In 
this region the ratio c^/v\ is equal to 1.4, 0.46 and 0.11 
in the three models with 5phot = 0.9, 1.5 and 2.4 kG, 
correspondingly. Thus, the wave source is located in the 
region dominated by the magnetic field in the last two 
models. In the second series, the source is located in the 
"penumbral" part of the model, at a distance of 12 Mm 
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from the axis, where the field is weaker. There the pa- 
rameter (?slv\ is equal to 8.2, 3.1 and 1.1 in the three 
models. In this case the waves are excited in the acous- 
tically dominated region. Placing the wave source at 
these two positions allows us to study the wave proper- 
ties for the two different regimes. Note, that through the 
whole paper we choose to discuss the results in terms 
of the parameter c|/v^ instead of plasma f5. The for- 
mer parameter is more relevant regarding wave propaga- 
tion. The ratio between these two parameters is equal 
(^I/'^a)//^ ~ 7/^ ~ ^/^ (fo^ ideal gas). The layers 
c|/v^ = 1 and = 1 are located very close in the atmo- 
sphere. The maximum height separation between them 
is about 30 km. The locations of the sources in both sets 
of simulations are marked by yellow dots in Fig. [TJ The 
duration of each simulation run is about 50 minutes. 

4.1. Non-magnetic case 

Fig. [3] shows the time evolution of the simulated hori- 
zontal and vertical velocities in the case of the quiet Sun 
model. In the figure, the velocities are scaled with the 
square root of the background density in order to make 
visible variations in the deep layers. The simulation re- 
veals that the single acoustic source excites a mixture 
of acoustic (p), gravity (g) and surface (/) modes, each 
propagating with its own speed. The p modes are the 
fastest and run beyond the boundaries of the simulation 
domain after about 30 min. The slower g and / modes 
are concentrated in the upper part of the simulation do- 
main. It can be seen that the wave field is symmetric 
with respect to the source location since the medium is 
isotropic for the acoustic-gravity waves. 

The power spectrum of the photospheric vertical ve- 
locity (k — v diagram) in the non-magnetic simulations 
is given in Fig. [H The diagram demonstrates the pres- 
ence of several ridges of acoustic (p-mode) oscillations. 
The darkest bottom ridge corresponds to / mode. The 
/-mode analytical dispersion relation is also plotted in 
the image and shows a reasonable agreement with the 
numerical results. The figure also reveals the presence 
of a weak ^-mode signal below the /-mode ridge. These 
modes appear because the background model has been 
made convectively stable. 

4.2. Waves in sunspot models 

Figures [5] and [6] show the results of the simulations in 
the sunspot model with 5phot = 1.5 kG and the source 
location at Xq = 6 Mm and 12 Mm, correspondingly. 
In the case of the simulations with magnetic field, we 
use for illustrations the longitudinal (parallel to the local 
field) and transversal (perpendicular to the field) velocity 
components. This provides an effective separation of the 
fast and slow magneto-acoustic waves in those regions of 
the simulation domain where the characteristic propaga- 
tion speeds are clearly distinct {va <^ cs or va ^ cs). 
Apart from a smaller propagation speed, the slow mode 
is mainly magnetic in the region where va <C c^, and has 
a displacement velocity that is perpendicular to the mag- 
netic field lines. This means that it contributes mainly 
to the transversal velocity. On the other hand, the fast 
mode is mainly acoustic, i.e., longitudinal, by its nature. 
As this can be seen in Fig. [5l the longitudinal compo- 
nent of the fast mode vanishes approximately along the 



field line crossing the source. The direction of propaga- 
tion of the slow magneto-acoustic waves (defined as the 
direction of their group velocity) is always parallel to the 
magnetic field direction, irrespective of their phase ve- 
locity. The direction of propagation of the fast magneto- 
acoustic waves is independent of the field direction, and 
the group velocity is parallel to the phase velocity. Given 
these properties, the fast and slow modes can be distin- 
guished in the time evolution pictures of the transversal 
and longitudinal velocity components, displayed in Fig- 
ures [5] and [6l 

In all simulations, the wave source is located in the re- 
gion where c'g/v\ stay (roughly) between 0.1 and 10, and 
the field is significantly inclined. In this region, both the 
pressure and magnetic restoring forces are of the same or- 
der of magnitude. Thus, the source generates a set of fast 
and slow magneto-acoustic waves. The fast modes are an 
analog of the acoustic p-modes in the non-magnetic sim- 
ulations (see Fig. [3]). However, as it will be shown later 
in the paper, their behavior is modified by the presence 
of the magnetic field. The slow modes did not appear in 
the non-magnetic simulation. 

Fig. [5] shows that the fast wave propagating in this 
complex magnetic field configuration causes variations 
in both, longitudinal and transversal, components of the 
velocity of the same order of magnitude. The initial wave 
front expands propagating downwards. There is no sym- 
metry in the wave field with respect to the source location 
(unlike the case of the non-magnetic simulation) since the 
wave propagation speed now depends on the direction. In 
the deep layers, the dominant velocity variations are due 
to the fast mode waves. In the upper photospheric layers 
we observe mostly variations in the longitudinal velocity 
(right panels in Fig. [5j). The transversal velocity varia- 
tions, caused by the fast magnetic waves {cg/v\ <C 1), 
are negligible in these layers. The wave velocity field 
developed in the simulations with the source located at 
Xq = 6 and 12 Mm from the sunspot axis are similar. 
The differences are quantitative, not qualitative in na- 
ture. 

Fig.[6]gives a more detailed view on the evolution of the 
transversal velocity component in the upper part of the 
simulation domain (note that the vertical axis has been 
expanded) and allows to understand better the complex 
physics of the wave propagation in the vicinity of the 
cs = '^A layer. An even better view on the wave dynam- 
ics in this region can be gained by watching the movie 
of the simulation, provided as on-line material to this 
paper. We observe in Fig. [6] the presence of the slow 
MHD mode generated directly by the source. This mode 
propagates with a visibly low speed downwards along the 
sunspot magnetic field lines. It is localized close to the 
horizontal position X = — 10 Mm. 

The fast waves generated by the source experience a 
significant refraction and refiection due to the strong in- 
crease of the Alfven speed (see Fig. [2]) in the vicinity 
of the cs = Va layer (snapshots corresponding to times 
between 8 and 14 minutes in Fig. [6]). 

In addition to the slow MHD waves generated directly 
by the source there is another wave type. It can be ob- 
served in Fig. [5] between 14 and 24 minutes of time as a 
perturbation with much smaller vertical wavelength com- 
pared to that of the fast modes in the sunspot interior. 
We can observe from Fig. [5]that the propagation speed of 
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Fig. 3. — Time series of snapshots of the horizontal (left panels) and vertical (right panels) velocities in the non-magnetic simulations 
with Xq = 6 Mm. The velocities are scaled with the color scale is the same in both panels. The negative heights correspond to 

sub-phot ospheric layers. 



tliis low-wavelengtli disturbance is comparable to that of 
the fast modes. Unlike the slow MHD waves, these waves 
propagate horizontally across the sunspot. Thus, these 
waves are different from either fast or slow MHD waves in 
several aspects: (1) they have a vertical wavelength that 
is much smaller than that of the fast waves; (2) they 
propagate with speeds similar to the fast waves (note, 
however, that they propagate in the region where the 



sound speed is slightly larger than the Alfven speed); (3) 
they propagate across the field, a behavior impossible for 
slow waves in this region; (4) the waves that posses sim- 
ilar characteristics (superficial horizontal propagation, 
smaller wavelength, etc) in non-magnetic simulations are 
the /-mode waves. From these arguments we suggest 
that we see a wave type with properties that are a mix- 
ture of magnetic and /-mode waves. We call them "sur- 



6 



Khomenko et al. 



p modes 










f mode 

g modes 





1 2 3 4 5 

\ [Mm-'] 

Fig. 4. — The k — v power spectrum of the photospheric vertical 
velocity in the non-magnetic standard solar model simulation run. 
Darker colors indicate larger power. The black line indicates the 
dispersion relation of the /-mode: uP' = gkx, where g is the value 
of solar gravity acceleration at the surface. Several p modes can be 
distinguished above the / mode. The g modes, located below the 
/ mode are also excited by the vertical force source. These modes 
are present because the solar model has been made convectively 
stable. 



face magneto-gravity" w aves. These wave s were men- 
tioned previously by, e.g. ISchwartz fc SteinI (|l975l ). who 
described their properties in an isothermal atmosphere. 
The variations produced by the surface magneto-gravity 
waves decrease rapidly with depth and disappear almost 
completely below —3 Mm. 

The power spectra of the photospheric vertical veloc- 
ity in the magnetic simulations for different 5phot and 
the source location at 12 Mm from the sunspot axis are 
compared in Fig. H Figure [3 shows that for the mag- 
netic models the k — u diagrams become progressively 
more and more different from the non-magnetic /c — di- 
agram (Fig. [4|), with increasing magnetic field strength. 
Even in the case of the weakest field, 5phot = 0.9 kG, 
the inclination of the p-mode ridges is larger than in the 
pure acoustic case. The inclination of the ridges increases 
with increasing magnetic field strength. The dispersion 
relation for the p-mode ridges depends on the reflection 
conditions for waves in the upper layers. In the absence 
of magnetic field, the waves with different frequencies 
reflect at different heights in the photosphere and be- 
low due to the rapid change of the cut-off frequency with 
height. When the magnetic field is introduced, the upper 
reflection condition changes even if the field is as weak 
as in the 5phot = 0.9 kG model. This happens because 
the propagation path of the high frequency fast waves 
is modified strongly by the rapid increase of the Alfven 
speed since these waves penetrate higher up in the at- 
mosphere. The reflection height for such waves becomes 
lower compared to the non-magnetic case. This naturally 
produces the change of the dispersion relation observed 
at the k — u diagrams in Fig. [71 

4.3. Time- distance diagrams and travel time estimates 



In order to quantify the effects of magnetic field on 
the travel time measurements of local helioseismology, we 
have constructed time-distance diagrams for the vertical 
velocity. Such diagrams were constructed for the photo- 
spheric velocity measured at the optical depth logrs = 
— 1.6 (see Sect. 3) in the sunspot models. Since in our 
simulations the waves are excited by a single source, 
there is no need to compute the velocity correlations 
between the different surface points. Thus, the time- 
distance wave propagation diagrams are constructed by 
plotting the velocity variations as a function of time for 
all horizontal points at constant logrs level. Fig. [8] shows 
the time-distance diagrams in the non-magnetic simula- 
tions and the simulations with progressively increasing 
magnetic field strength. 

In the time-distance diagram for the non-magnetic case 
(Fig. [HI a), the two lowest system of ridges are due to the 
p modes. The green dashed line in Fig. [8^ gives a the- 
oretical dispersion relation for these modes, calculated 
for the standard solar model S. Above them, the third 
system of ridges comes from the / mode waves. It can 
be seen from the inclination of the ridges that the phase 
and group speeds of this mode are different from those of 
the p modes, confirming that they belong to the /-mode. 
The theoretical time-distance curve for the /-mode (blue 
dashed line in Fig. [8tt) is in a good agreement with the 
numerical results. The system of ridges at the top part 
of the diagram above the source corresponds to the g 
modes. These modes propagate with much lower speed 
than the p and / modes, in agreement with their appear- 
ance at the k — u diagram in Fig. [H 

The time-distance diagrams in the simulations with 
magnetic field have different properties (Fig. [51 b, c, and 
d). p- mode-like ridges are present in all the cases. As in 
the non-magnetic case, we have added theoretical curves 
for p-modes and the /-mode in the absence of the mag- 
netic field (green and blue dashed curves, in Fig. [Hb, c, 
and d). These curves help identifying the fast modes 
as the analog of the non-magnetic p-modes propagat- 
ing with a speed modified by the magnetic field. The 
/- and ^-mode ridges do not appear, in agreement with 
the k — V diagram in Fig. [71 Instead, weak traces of 
the slow waves are present in the top part of the dia- 
gram. This mode is not symmetric with respect to the 
source location and has only surface signatures on the 
side of the source where the magnetic field is larger. At 
the logTs = —1.6, the slow mode is mostly prominent in 
the sunspot simulation with weak field 5phot = 0.9 kG. 
It becomes weaker but still can be distinguished in the 
models with stronger field. The origin of this mode is 
different from the slow modes generated by the source in 
the sunspot interior (Fig. [6]) and needs further study. 

At first glance, the p-mode ridges on the time-distance 
diagrams seem not to be affected much by the magnetic 
field. A more detailed investigation shows that the p- 
mode ridge, at the bottom, becomes thicker, i.e. the 
relative contribution of the different frequencies in the 
wave packet changes with magnetic field. We have per- 
formed a fit to the time-distance diagrams in order to 
find the phase and group travel times of waves in each 
case. At each horizontal position starting from 3 Mm 
from the source we fit the time velocity curve with a 
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Fig. 5. — Time series of snapshots of the transversal (left panels) and longitudinal (right panels) velocities in the simulations with Bphot 
=1.5 kG and Xq = 6 Mm. The velocities are scaled with the factor of ^/pQ. The color scales are the same in all panels. The red and 
green lines are contours of constant c^s/^\^ ^^^^ corresponds to va = cs, and the green line corresponds to c^/i;^ =0.1. The black 

inclined lines are the magnetic field lines. The negative heights correspond to sub-photospheric layers. Note the presence of slow waves, in 
addition to the modified p-modes (fast waves), indicated by arrows in the time = 15 min panel. 



Gabor's wavelet: 

G(r) = Acos[u;o{r — Tp)] exp 



4 



(8) 



where A is the amplitude, ujq is the central frequency, Tp 
and Tg are the phase and group travel times respectively. 



and SiJ is the bandwidth fl Kosovichev fc Duvalll 119971 ). 
The fit has four free parameters: cc^o, r^, Tg^ and 6uj. The 
red curves plotted in the diagrams of Fig. [S] are the phase 
travel times Tp obtained after the fit. 
Fig. [9] shows the differences in the phase travel-times 
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Fig. 6. — Time series of snapshots of the transversal velocity in 
the simulations with Bp^ot =1-5 kG and Xq = 12 Mm. Only the 
top part of the simulation box is shown and the vertical axis has 
been expanded for better presentation. The velocities are scaled 
with the factor of ^/po. Each panel is normalized to its own maxi- 
mum. 



between the magnetic and non-magnetic simulations. 
Negative values mean that waves in the magnetic simula- 
tions propagate faster. The travel- time difference curves 
have a complex dependence on the distance. The dif- 
ference between the magnetic and non-magnetic cases 
is larger at locations closer to the sunspot axis, where 
the magnetic field is stronger, and getting smaller at 
15 Mm far from the axis. This happens because the 
waves with larger skip distances travel deeper inside the 
sunspot and are less affected by the magnetic surface 
effects. The travel time difference depends on the mag- 
netic field strength of the model. The model with the 
smallest field 5phot = 0-9 kG gives a maximum differ- 
ence in the travel times of about 25 seconds, while in 
the model with the largest field strength 5phot = 2.4 



kG, this difference becomes about 35 seconds. The re- 
sults presented in Figs. [8] and [9] are calculated for the 
photospheric vertical velocities taken at constant logrs 
level. These results do not change much if we take ve- 
locities at a constant geometrical height Z = 100 km, 
instead of constant logr5=— 1.6. Only a minor additional 
phase shift is introduced in the travel time variations in 
Fig. [9] making them about 5 seconds smaller close to the 
sunspot axis. 

The travel-time differences discussed above were 
calculated for the central frequency of the wave 
packet. However, there has been a discussion in 
the recent literature regarding the frequency depen- 
dence of the travel time me asurements (Couvid at et al.l 
120061: 'Couvidat & Raiagurul 120071 : iRajaguru "2008; 
iBraun & Birch 2008). It has been obtained that, if the 
frequency filtering is applied to the data, in addition to 
the phase speed filtering, the resulting travel-time dif- 
ferences show a strong frequency dependence. High fre- 
quenc y waves have, in general, larger travel-time differ- 
ences (jRajagurull2008[ ) . In order to investigate this effect 
in our simulations, we filtered the simulated photospheric 
vertical velocities in the five frequency ranges: 3.3, 3.5, 
4.0, 4.5 and 4.8 mHz. We applied a Gaussian filter with 
a FWHM of 1 mHz. The resulting travel-time differences 
are shown in Fig. [TOl for all three model sunspots. 

We find that the travel-time differences for each in- 
dividual frequency behave similarly to the travel-time 
differences calculated without filtering. The maximum 
difference is achieved for distances closer to the sunspot 
axis (±5 Mm), decreasing for the waves with large skip 
distances. A strong frequency dependence is evident 
from Fig. [TOl Waves with higher frequencies have larger 
travel- time differences. It means that the propagation 
speed of waves increases with frequency. The magni- 
tude of the frequency dependence that we obtain from 
the simulations is in agreement with the observations of 
iRaiagurul (|2008[ ). In that paper travel-time differences 
change from, about, —15 to —30 seconds between 3 and 
5 mHz, similar to our results in Fig. [TOl 

Interestingly, the travel-time differences get slightly 
positive for the distances within first 5 — 6 Mm from 
the source location and for lower frequencies between 3.3 
and 4.0 mHz. This effect needs further investigation. 

5. EIKONAL APPROXIMATION SOLUTIONS 

In order to understand more deeply the numerical re- 
sults, we present here an approximate solution of the 
governing MHD equations for small perturbations in 
the eikonal approximation (se e recent applications in 
e.g. [M cLaughlin & Hood 2006; i McLaughhn et al.ll2008l : 
Khomenko fc Colla dos 2006). The eikonal approach as- 
sumes that the wavelength of the perturbation is much 
smaller than the characteristic scale of the variations of 
the background model (^Gou gh*l2007l). In the zero-order 
eikonal approximation, we neglect the variation of the 
wave amplitude and consider only the variation of its 
phase, i.e., we assume that all the variables in ([l]-[4|) de- 
pend on X, z and time as /7 = a • e^^^^^'^^ • e"*'^^ (where 
a is constant). We ma ke use of the recent developments 
of lcidWl ([2006h and Mor adi CaTh^ (2008) that allow to 
incorporate the effects of the acoustic cut-off frequency 
in the eikonal solution for MHD waves. In this approxi- 
mation after some mathematical simplifications, the fol- 
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Fig. 7. — The k — u power spectra of the photospheric vertical velocity in the simulation runs for the source position Xq= 12 Mm and 
5phot=0-9 (Q'), 1-5 (b) and 2.4 (c) kG. Darker colors indicate larger power. Note the change of the power spectra with the increasing field 
strength. 



lowing first-order non-linear partial differential equation 
is obtained: 



(9) 



2/ 2 



2 2\ 

csa ) 



where P = and <7 = ff • Parameters p and q are equiv- 
alent to the horizontal and vertical wave numbers, re- 
spectively. Parameter N is the Brunt- Vaisala frequency, 
N'^ = g/H — g'^/c^. In the following calculations, we use 
an expression for the cut-off frequency for an isother- 
mal atmosphere ujc = cs/2H. Equation (p!Q|) is solved 
by Charpit's method of characteristics by transforming 
it into the following system of ordinary differential equa- 
tions: 



dp 

ds 
dq 

ds 
dx 

ds 
dz 
ds 



dF 
dx 

dz 
dF 

dp 
dF 



(10) 

(11) 
(12) 

(13) 



Here, variable s is the distance along the characteris- 
tic wave propagation path. These four ODEs are solved 
numerically using a fourth-order Runge-Kutta method. 
We investigate only fast-mode solutions of the equations. 
Parameters c^, va^ N and Uc are taken from the sunspot 
models. The solution of Eqs. (p!QUT3|) gives us x{s)^ ^{^)^ 
p{s) and q{s) along the wave path s. The lines x{z) 
give the trajectory of the group velocity of the solu- 
tio n. The phase alon^ th e wave path can be calculated 
as (|Moradi fc Callvll2008D : 



kdf — cot 



(14) 



where r = (x, z) and k = (p, q). The phase function S{f) 
divided by the wave frequency co gives us an estimate of 
the phase travel time. 

First, we study the paths through the sunspot models 
of fast magneto- acoustic waves of a given frequency co. 
Waves are launched from the same lower turning point 



at some depth below the surface. Fig. [TT] gives the exam- 
ples of solutions x{z) for the fast wave paths for the fre- 
quencies 6.6 (top), 4.5 (middle) and 3.0 (bottom) mHz. 
As initial condition for the solution of the equations (p^Ql - 
fT3|) we take in all the cases x = —40 Mm, z = —3.5 
Mm and q = 0. The initial value of p is found from the 
approximated solution of the dispersion relation (p!Q|) at 
q = 0. 

At each frequency, we calculate the wave paths for 
two cases. In the first case (red lines in Fig. [TT]) we set 
the Alfven speed va to zero {i.e. eliminate the magnetic 
field), but keep the sound speed structure and frequen- 
cies TV and cOc from the sunspot model. In the second 
case, we calculate the wave path in the sunspot model 
taking into account magnetic field terms (white lines in 
Fig. [TT]) . As can be observed in Fig.[TTl the paths of the 
waves in the magnetic and non-magnetic cases are sig- 
nificantly different. Black dots on the curves mark the 
positions separated 1 minute in time in order to outline 
the quantitative differences in phase times spent by the 
waves on their paths. 

The frequency of the wave at the top panel of Fig. (TTj 
(6.6 mHz) is larger than the atmospheric cut-off. Thus, 
in the non-magnetic case, this wave (red line) propagates 
to the upper atmospheric layers without reflection and 
goes outside the simulation domain. This does not hap- 
pen in the magnetic case. The green horizontal line in 
Fig.[TTlmarks the position of the layer where cs = va- As 
the fast magneto-acoustic wave reaches this layer from 
below, its propagation properties become progressively 
more affected by the magnetic field. The fast magneto- 
acoustic wave in the region va > cs experiences the re- 
fraction due to the rapid increase of the va with height (a 
clear example of s uch behavior can b e seen in the simula- 
tions presented bv lKhomenko fc Co llados 2006). Finally 
the fast magneto-acoustic wave reflects back to the deep 
layers at some height above the cs = va- In the deep 
layers, this wave experiences again the refraction and re- 
flection, but now due to the rapid increase of sound speed 
Cs- Thus, in this model, the high-frequency fast waves 
in the magnetic case are trapped below the photosphere 
due to magnetic effects. 

If we reduce the wave frequency to 4.5 mHz (middle 
panel in Fig. [TT]) . then in both the magnetic and non- 
magnetic cases the wave paths have an upper reflection 
point. The height where the wave frequency is equal to 
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Fig. 8. — The time-distance diagrams of the photospheric vertical 
velocity for the simulations with (a) no magnetic field; (b) sunspot 
model with Sphot = 0-9 kG; (c) sunspot model with Sphot = 1-5 
kG; (d) sunspot model with Bp^ot = 2.4 kG. The source is located 
at Xq = 12 Mm on the left of the sunspot axis. The red curves are 
the phase travel times calculated using the Gabor wavelet fit to the 
time-distance curves (see the text). The green dashed curves are 
theoretical time-distance curves for p-modes and the standard solar 
model S in the absence of the magnetic field. The blue dashed lines 
are the theoretical time-distance curves for the / mode calculated 
as X = (l/2)(g/u;)t for the central frequency u; = 27r/300 of the 
wave source. 



the atmospheric cut-off frequency is marked by the yehow 
hues in the figure. Since the sunspot model atmosphere 
has a temperature minimum, there is a region in between 
the two yehow hues, where the wave of frequency 4.5 
mHz is evanescent. In the non- magnetic case, the wave is 
reflected sharply from this region. In the magnetic case, 
the reflection process is more complex. One can observe 
that the relative positions of the layer where cs = va and 
the cut-off {uJc) layer change with the horizontal distance. 
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Fig. 9. — The phase travel-time difference calculated between 
the travel times measured in simulations without magnetic field 
and the simulations with different model sunspots (indicated in the 
figure) as a function of horizontal distance. The source location at 
Xq = 12 Mm left to the sunspot axis is marked by yellow circle. 



Far from the sunspot axis, where the field is relatively 
weak, the cut-off height lies below the cs — va layer. 
Thus, the waves are reflected before magnetic effects on 
the wave path take place. Note that the trajectories of 
the waves in the magnetic and non-magnetic cases are 
very similar in this region. As the waves approximate the 
region close to the sunspot center (where the magnetic 
field is large), the position of the cs = va layer lowers 
and the fast wave path becomes affected by the magnetic 
field. It can be seen from the figure that the red and 
white curves start to separate significantly in the region 
close to the sunspot axis. The waves in the magnetic case 
propagate with larger speeds and overtake the waves in 
the non-magnetic case. Note that by introducing the 
magnetic field, we have changed the height of the upper 
turning point of the waves. This may explain the shift 
of the positions of the ridges in the magnetic simulations 
in Fig. [71 compared to the non-magnetic simulations in 
Fig.H 

Finally, the bottom panel of Fig. [11] gives the propa- 
gation paths in the low- frequency case of 3 mHz. In this 
case, the position of the cut-off height always lies below 
the height of cs = va layer. It means that (at least in 
our sunspot model atmosphere) low frequency waves are 
reflected before their path is significantly affected by the 
rapid increase of the Alfven speed. However, it is im- 
portant to note that even in this case, the trajectories of 
the waves in the magnetic and non-magnetic models are 
different. The propagation speed of the fast magneto- 
acoustic wave depends on the magnetic field and always 
greater than the sound speed. 

6. QUANTITATIVE ESTIMATES OF THE TRAVEL TIMES 
AND SKIP DISTANCES 

In order to quantify the effects of the magnetic field on 
the wave travel times in the eikonal approximation, we 
calculated the differences in travel times for the waves in 
the magnetic and non-magnetic cases. This calculation is 
done for waves frequencies between 3.3 and 4.8 mHz. The 
resulting travel time differences are presented in Fig. [12] 
as a function of the wave skip di stance. This calculatio n 
is similar to the one reported bv lMoradi fc Callvl (|2QQ8[ ). 
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Fig. 10. — The phase travel time difference between the magnetic 
and non-magnetic cases as a function of the horizontal distance, 
calculated for different frequency filtering of the simulations with 
the source located at 12 Mm. The different colors correspond to 
different wave frequencies. Calculations done for the model sunspot 
with Bphot = 0.9 kG (top), 1.5 kG (middle) and 2.4 kG (top). 



The curves shown in Fig. [12] were calculated in the fol- 
lowing way. In each case, the waves are launched from 
their lower turning point, thus g = 0, and p is obtained 
from the dispersion relation (p!Q|) . The lower turning 
point is located at annuli distance xa = 12 Mm from 
the sunspot axis, in order to make the eikonal calcula- 
tions directly comparable to the travel time estimates 
from simulations, as presented in Fig. [TOl The vertical 
location of the lower turning point za varies in a range 
from —1.5 to —10 Mm below the photosphere. By vary- 
ing Za and frequency cj, we obtain different curves in 
Fig. [121 The deeper is za, the larger is the value of the 
skip distance in the horizontal axis of Fig. [121 Note that 
to facilitate the comparison between the numerical sim- 



ulations and the eikonal approach, we present the skip 
distances in Fig. [121 with respect to the position to the 
sunspot axis. 

The frequency behavior calculated in the eikonal ap- 
proximation is in good agreement with the simulations. 
It proves the reliability of the eikonal approximation 
for the full MHD approach. As in Fig. [TQl the travel 
time differences are negative except for positions close 
to the source location, meaning that the waves propa- 
gate faster in the magnetic case. There is a strong fre- 
quency dependence of the travel time differences. Fig. [121 
shows that depending on the sunspot model and on the 
skip distance, the maximum values change from, about, 
—50 sec for the 4.8 mHz wave to —10 sec for the 3.3 
mHz wave. The agreement between simulations and 
the eikonal solution is worse in the 5phot = 2.4 kG 
case, where the latter gives slightly larger values of the 
travel time differences and a different dependence. In 
all the cases the radius dependence obtained from the 
simulations using the Gabor wavelet fit is more compli- 
cated than the skip-distance dependence obtained in the 
eikonal approximation. Note that both simulations and 
the eikonal approximation give values of the travel time 
variations that fit reasonably well into the range of the 
observed values in solar active regions below sunspots 
(see, e.g. Kosovichev et al. 2000; Zhao & KosovichevJ 
20031: ICouvidat et al 1 120061: ICouvidat Raiagurull2007l : 
Braun fc Birchll2008l ). However, the frequency depen- 
dence in the eikonal approximation is slightly larger than 
that obtained from simulations and observations. 

7. DISCUSSION AND CONCLUSIONS 

In this paper, we present an investigation of the influ- 
ence of magnetic field on helioseismic wave propagation 
in active regions below sunspots. In order to better un- 
derstand the physics of the waves in complex magnetic 
field configurations, we perform a single source experi- 
ment and study the propagation properties of the dif- 
ferent wave modes excited by such a source, as well as 
the wave travel times. We have carried out our calcu- 
lations in a series of sunspot models with different field 
strength and for two source positions located in magnet- 
ically and acoustically dominated regions. The results of 
our 2D simulations are compared to an asymptotic solu- 
tion using the eikonal approximation. The results of both 
methods, numerical simulations and eikonal approxima- 
tion, are found to be in a good qualitative agreement. 
Our approach has allowed us to investigate the effects 
of the magnetic field strength in sunspots on the travel- 
time variations obtained there from observations using 
the time-distance helioseismology techniques. In addi- 
tion, the frequency dependence of such measurements 
has been studied. 

The following conclusions summarize our study: 

• The wave source located immediately below the 
surface in a non-magnetic standard solar model ex- 
cites a mixture of p, / and g mode waves. When 
the magnetic field is present, the same source ex- 
cites a mixture of fast and slow MHD waves. The 
fast modes represent an analog of the acoustic p 
modes, modified by the magnetic field. 

• In addition to slow MHD waves excited directly 
by the source, there is another wave type with low 
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Fig. 11. — The wave paths of the fast mode waves launched from the same lower turning point below the photosphere propagating 
through the model sunspot with Bp^ot = 2.4 kG for different frequencies: 6.6 mHz (top); 4.5 mHz (middle) and 3 mHz (bottom). The 
white curves are the wave paths when the magnetic field is taken into account. The red curves are the wave paths in the sunspot models 
with the same thermal properties, but setting B = 0. The black dots on the trajectory are separated from each other by 1 minute in time. 
The green curve marks the position of the va = cs layer. The yellow curves mark the heights where the frequency of the wave is equal to 
the local acoustic cut-off frequency. The background grey image is acoustic speed cs- 



vertical wavelength, propagating with speed close 
to that of the fast waves horizontally across the 
sunspot. We call them surface magneto-gravity 
waves. Several arguments listed in Sect. 4 indi- 
cate that these waves are different from either fast 
or slow waves. In our simulations these waves prop- 
agate below the visible surface. This makes their 
detection complicated in spectral observations. 

• Slow MHD wave modes form an additional group 
of ridges in the time-distance diagram with much 
lower propagation speed. Their properties depend 
on the magnetic field strength of the model, as well 
as on the height in the atmosphere where the ve- 
locity is measured. The slow mode ridges on the 
time-distance diagrams are well isolated and do not 
affect the travel time measurements in magnetic re- 
gions. 

• The helioseismic waves below sunspots are speed 
up by the magnetic field. The travel times of these 
waves are shorter by, about 20-40 seconds com- 



pared to the waves in the quiet Sun. Changing the 
photospheric magnetic field strength from 0.9 to 
2.4 kG produces a variation of the travel times by 
about 10-15 seconds. 

• Magnetic field produces a strong frequency depen- 
dence of the wave travel times. High frequency 
waves travel faster than low frequency waves. This 
happens because the propagation path of high fre- 
quency fast waves is modified strongly by the rapid 
increase of the Alfven speed in the top layers, since 
these waves penetrate higher up in the atmosphere. 
Low frequency waves are refiected below the pho- 
tosphere due to the acoustic cut-off frequency and 
do not reach the heights where va > cs and the 
magnetic fields produces strong effects. The travel 
time differences change from —20 to —35 seconds 
between 3 and 5 mHz, in agreement with observa- 
tions (e.g. iRa jagur ^ l2008h . 

The analysis presented here is limited by the 2D ap- 
'oximation. As a next step in our work we plan to relax 
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this assumption and to analyze 3D simulations, with a 
wave excitation by either single or multiple sources. This 
will allow us to build a more realistic picture and to com- 
plete our understanding of the interaction of helioseismic 
waves with magnetic fields of sunspots. 
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Fig. 12. — The travel time difference between magnetic and 
non-magnetic cases as a function of the skip distance, calculated in 
the eikonal approximation for the annuli size of 12 Mm. The lines 
of different colors correspond to different wave frequencies. The 
calculations are done for the sunspot models with Bp^ot =0.9 kG 
(top), 1.5 kG (middle) and 2.4 kG (top). 
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